Loading required packages
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(tidyverse)
library(sf)
library(tmap)
library(PerformanceAnalytics) #scatter plot matrix
library(PerformanceAnalytics) #scatter plot matrix
# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=5, scipen=999)
Import the shapefile
sh<- st_read("../../Data/city.shp")
Reading layer `City' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\book-2020-spatial-analysis-methods-and-practice\Data\City.shp' using driver `ESRI Shapefile'
Simple feature collection with 90 features and 15 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 472360 ymin: 4199900 xmax: 481460 ymax: 4209300
projected CRS: GGRS87 / Greek Grid
140128 bytes
s.sp <- as(sh, "Spatial")
class(s.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
[1] "sf" "data.frame"
Inspect the table within the shapefile
A quick map

Creating the breaks for the custom breaks map
breaks<- c(-Inf, 15000,20000,25000, Inf)
creating the custom breaks map
tm_shape(sh) +
tm_fill("Income", title= "Income",breaks=breaks,palette="Greens")+
tm_borders()

Creating a histogram
sh %>%
ggplot(aes(x=Income))+
geom_histogram(bins=10, color="black", fill="white")

creating a box plot
geom_jitter(size = 3, alpha = 0.3, width = 0.1) is used so that the points are not all crammed in the same line
set.seed(2019)
sh %>%
ggplot(aes(x=Income, y=Municipali))+
geom_boxplot()+
geom_jitter(size = 3, alpha = 0.3, width = 0.1)

Creating a new variable standardized income
sh$IncZScore<- as.vector(scale(sh$Income, center = TRUE, scale = TRUE))
creating the new breaks
breaks_Z<- c(-Inf, -2.5,-1, 1, 2.5, Inf)
creating the custom breaks map
#creating a vector with the colors for the bins in the map
col<- c('#ffffb2','#fecc5c','#fd8d3c','#f03b20','#bd0026')
tm_shape(sh) +
tm_fill("IncZScore", title= "Standardized Income", breaks=breaks_Z, palette = col)+
tm_borders()

Bivariate Analysis: Analyzing Expenditures by Educational Attainment
A scatter plot and the OLS regression estimates
sh %>%
ggplot(aes(x=University, y=Expenses))+
geom_point()+
geom_smooth(method = "lm")

lm.sh<- lm(Expenses ~ University, data=sh)
summary( lm.sh)
Call:
lm(formula = Expenses ~ University, data = sh)
Residuals:
Min 1Q Median 3Q Max
-132.33 -29.14 -14.76 9.47 293.99
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.69 19.17 0.82 0.42
University 8.89 0.97 9.16 0.000000000000019 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 64.4 on 88 degrees of freedom
Multiple R-squared: 0.488, Adjusted R-squared: 0.482
F-statistic: 83.9 on 1 and 88 DF, p-value: 0.0000000000000193
Creating a scatter plot matrix
data<- data.frame( sh$Expenses, sh$University, sh$SecondaryE)
data
chart.Correlation(data, histogram = TRUE, method = "pearson")

A better option
#pairs.panels(data,
# smooth = TRUE, # If TRUE, draws loess smooths
# scale = FALSE, # If TRUE, scales the correlation text font
# density = TRUE, # If TRUE, adds density plots and histograms
# ellipses = TRUE, # If TRUE, draws ellipses
# method = "pearson", # Correlation method (also "spearman" or "kendall")
# pch = 21, # pch symbol
# lm = TRUE, # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
# cor = TRUE, # If TRUE, reports correlations
# jiggle = FALSE, # If TRUE, data points are jittered
# factor = 2, # Jittering factor
# hist.col = 4, # Histograms color
# stars = TRUE, # If TRUE, adds significance level with stars
# ci = TRUE) # If TRUE, adds confidence intervals
That’s all!

LS0tDQp0aXRsZTogIkV4cGxvcmF0b3J5IFNwYXRpYWwgRGF0YSBBbmFseXNpcw0KVG9vbHMgYW5kIFN0YXRpc3RpY3MiDQphdXRob3I6ICJGZWxpcGUgU2FudG9zLU1hcnF1ZXogJiBTZXcgU29vayBZYW4iDQpkYXRlOiAiMTAvMTIvMjAyMCINCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazoNCiAgICBjb2RlX2ZvbGRpbmc6IHNob3cNCiAgICBoaWdobGlnaHQ6IG1vbm9jaHJvbWUNCiAgICBudW1iZXJfc2VjdGlvbnM6IHllcw0KICAgIHRoZW1lOiBjb3Ntbw0KICAgIHRvYzogeWVzDQogICAgdG9jX2RlcHRoOiA0DQogICAgdG9jX2Zsb2F0Og0KICAgICAgY29sbGFwc2VkOiBubw0KICAgICAgc21vb3RoX3Njcm9sbDogbm8NCiAgaHRtbF9kb2N1bWVudDoNCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQogICAgZGZfcHJpbnQ6IHBhZ2VkDQogICAgdG9jOiB0cnVlDQogICAgdG9jX2Zsb2F0Og0KICAgICAgY29sbGFwc2VkOiBmYWxzZQ0KICAgICAgc21vb3RoX3Njcm9sbDogZmFsc2UNCiAgICB0b2NfZGVwdGg6IDQNCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUNCiAgICBjb2RlX2ZvbGRpbmc6ICJzaG93Ig0KICAgIHRoZW1lOiAiY29zbW8iDQogICAgaGlnaGxpZ2h0OiAibW9ub2Nocm9tZSINCiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0DQogIHdvcmRfZG9jdW1lbnQ6IGRlZmF1bHQNCmJpYmxpb2dyYXBoeTogYmlibGlvLmJpYg0KLS0tDQoNCjxzdHlsZT4NCmgxLnRpdGxlIHtmb250LXNpemU6IDE4cHQ7IGNvbG9yOiBEYXJrQmx1ZTt9IA0KYm9keSwgaDEsIGgyLCBoMywgaDQge2ZvbnQtZmFtaWx5OiAiUGFsYXRpbm8iLCBzZXJpZjt9DQpib2R5IHtmb250LXNpemU6IDEycHQ7fQ0KLyogSGVhZGVycyAqLw0KaDEsaDIsaDMsaDQsaDUsaDZ7Zm9udC1zaXplOiAxNHB0OyBjb2xvcjogIzAwMDA4Qjt9DQpib2R5IHtjb2xvcjogIzMzMzMzMzt9DQphLCBhOmhvdmVyIHtjb2xvcjogIzhCM0E2Mjt9DQpwcmUge2ZvbnQtc2l6ZTogMTJweDt9DQo8L3N0eWxlPg0KDQojIExvYWRpbmcgcmVxdWlyZWQgcGFja2FnZXMNCg0KYGBge3Igc2V0dXB9DQprbml0cjo6b3B0c19jaHVuayRzZXQod2FybmluZyA9IEZBTFNFLCBtZXNzYWdlID0gRkFMU0UsIGVjaG8gPSBGQUxTRSkNCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShzZikNCmxpYnJhcnkodG1hcCkNCmxpYnJhcnkoUGVyZm9ybWFuY2VBbmFseXRpY3MpICNzY2F0dGVyIHBsb3QgbWF0cml4DQpsaWJyYXJ5KFBlcmZvcm1hbmNlQW5hbHl0aWNzKSAjc2NhdHRlciBwbG90IG1hdHJpeA0KIA0KIyBDaGFuZ2UgdGhlIHByZXNlbnRhdGlvbiBvZiBkZWNpbWFsIG51bWJlcnMgdG8gNCBhbmQgYXZvaWQgc2NpZW50aWZpYyBub3RhdGlvbg0Kb3B0aW9ucyhwcm9tcHQ9IlI+ICIsIGRpZ2l0cz01LCBzY2lwZW49OTk5KQ0KYGBgDQoNCg0KIyBJbXBvcnQgdGhlIHNoYXBlZmlsZQ0KDQoNCmBgYHtyfQ0Kc2g8LSBzdF9yZWFkKCIuLi8uLi9EYXRhL2NpdHkuc2hwIikNCm9iamVjdC5zaXplKHNoKQ0Kcy5zcCA8LSBhcyhzaCwgIlNwYXRpYWwiKQ0KY2xhc3Mocy5zcCkNCmNsYXNzKHNoKQ0KYGBgDQoNCiMjIEluc3BlY3QgdGhlIHRhYmxlIHdpdGhpbiB0aGUgc2hhcGVmaWxlDQoNCmBgYHtyfQ0Kc3RfZHJvcF9nZW9tZXRyeShzaCkNCmBgYA0KDQojIyBBIHF1aWNrIG1hcA0KDQpgYGB7cn0NCnF0bShzaCkNCmBgYA0KIyBDcmVhdGluZyB0aGUgYnJlYWtzIGZvciB0aGUgY3VzdG9tIGJyZWFrcyBtYXAgDQoNCmBgYHtyfQ0KYnJlYWtzPC0gYygtSW5mLCAxNTAwMCwyMDAwMCwyNTAwMCwgSW5mKQ0KYGBgDQoNCiMjIGNyZWF0aW5nIHRoZSBjdXN0b20gYnJlYWtzIG1hcA0KDQpgYGB7cn0NCnRtX3NoYXBlKHNoKSArDQogIHRtX2ZpbGwoIkluY29tZSIsIHRpdGxlPSAiSW5jb21lIixicmVha3M9YnJlYWtzLHBhbGV0dGU9IkdyZWVucyIpKw0KICB0bV9ib3JkZXJzKCkNCmBgYA0KDQoNCiMgQ3JlYXRpbmcgYSBoaXN0b2dyYW0NCg0KYGBge3J9DQpzaCAlPiUgDQogIGdncGxvdChhZXMoeD1JbmNvbWUpKSsNCiAgZ2VvbV9oaXN0b2dyYW0oYmlucz0xMCwgY29sb3I9ImJsYWNrIiwgZmlsbD0id2hpdGUiKQ0KYGBgDQoNCg0KIyBjcmVhdGluZyBhIGJveCBwbG90DQoNCmdlb21faml0dGVyKHNpemUgPSAzLCBhbHBoYSA9IDAuMywgd2lkdGggPSAwLjEpIGlzIHVzZWQgc28gdGhhdCB0aGUgcG9pbnRzIGFyZSBub3QgYWxsIGNyYW1tZWQgaW4gdGhlDQpzYW1lIGxpbmUNCg0KYGBge3J9DQpzZXQuc2VlZCgyMDE5KQ0Kc2ggJT4lIA0KICBnZ3Bsb3QoYWVzKHg9SW5jb21lLCB5PU11bmljaXBhbGkpKSsNCiAgZ2VvbV9ib3hwbG90KCkrDQogIGdlb21faml0dGVyKHNpemUgPSAzLCBhbHBoYSA9IDAuMywgd2lkdGggPSAwLjEpDQpgYGANCg0KIyBDcmVhdGluZyBhIG5ldyB2YXJpYWJsZSBzdGFuZGFyZGl6ZWQgaW5jb21lDQoNCmBgYHtyfQ0Kc2gkSW5jWlNjb3JlPC0gYXMudmVjdG9yKHNjYWxlKHNoJEluY29tZSwgY2VudGVyID0gVFJVRSwgc2NhbGUgPSBUUlVFKSkNCmBgYA0KDQojIyBjcmVhdGluZyB0aGUgbmV3IGJyZWFrcw0KDQpgYGB7cn0NCmJyZWFrc19aPC0gYygtSW5mLCAtMi41LC0xLCAxLCAyLjUsIEluZikNCmBgYA0KDQojIyBjcmVhdGluZyB0aGUgY3VzdG9tIGJyZWFrcyBtYXANCg0KYGBge3J9DQojY3JlYXRpbmcgYSB2ZWN0b3Igd2l0aCB0aGUgY29sb3JzIGZvciB0aGUgYmlucyBpbiB0aGUgbWFwDQpjb2w8LSBjKCcjZmZmZmIyJywnI2ZlY2M1YycsJyNmZDhkM2MnLCcjZjAzYjIwJywnI2JkMDAyNicpDQoNCnRtX3NoYXBlKHNoKSArDQogIHRtX2ZpbGwoIkluY1pTY29yZSIsIHRpdGxlPSAiU3RhbmRhcmRpemVkIEluY29tZSIsIGJyZWFrcz1icmVha3NfWiwgcGFsZXR0ZSA9IGNvbCkrDQogIHRtX2JvcmRlcnMoKQ0KYGBgDQoNCiMgQml2YXJpYXRlIEFuYWx5c2lzOiBBbmFseXppbmcgRXhwZW5kaXR1cmVzIGJ5IEVkdWNhdGlvbmFsIEF0dGFpbm1lbnQgDQoNCkEgc2NhdHRlciBwbG90IGFuZCB0aGUgT0xTIHJlZ3Jlc3Npb24gZXN0aW1hdGVzICANCg0KYGBge3J9DQpzaCAlPiUgDQogIGdncGxvdChhZXMoeD1Vbml2ZXJzaXR5LCB5PUV4cGVuc2VzKSkrDQogIGdlb21fcG9pbnQoKSsNCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikNCg0KbG0uc2g8LSAgIGxtKEV4cGVuc2VzIH4gVW5pdmVyc2l0eSwgZGF0YT1zaCkNCnN1bW1hcnkoIGxtLnNoKQ0KYGBgDQoNCiMgQ3JlYXRpbmcgYSBzY2F0dGVyIHBsb3QgbWF0cml4DQoNCmBgYHtyfQ0KZGF0YTwtIGRhdGEuZnJhbWUoIHNoJEV4cGVuc2VzLCBzaCRVbml2ZXJzaXR5LCBzaCRTZWNvbmRhcnlFKQ0KZGF0YQ0KY2hhcnQuQ29ycmVsYXRpb24oZGF0YSwgaGlzdG9ncmFtID0gVFJVRSwgbWV0aG9kID0gInBlYXJzb24iKQ0KYGBgDQoNCiMjIEEgYmV0dGVyIG9wdGlvbiANCg0KYGBge3J9DQojcGFpcnMucGFuZWxzKGRhdGEsDQogIyAgICAgICAgICAgIHNtb290aCA9IFRSVUUsICAgICAgIyBJZiBUUlVFLCBkcmF3cyBsb2VzcyBzbW9vdGhzDQogIyAgICAgICAgICAgIHNjYWxlID0gRkFMU0UsICAgICAgIyBJZiBUUlVFLCBzY2FsZXMgdGhlIGNvcnJlbGF0aW9uIHRleHQgZm9udA0KICMgICAgICAgICAgICBkZW5zaXR5ID0gVFJVRSwgICAgICMgSWYgVFJVRSwgYWRkcyBkZW5zaXR5IHBsb3RzIGFuZCBoaXN0b2dyYW1zDQogICMgICAgICAgICAgIGVsbGlwc2VzID0gVFJVRSwgICAgIyBJZiBUUlVFLCBkcmF3cyBlbGxpcHNlcw0KICAjICAgICAgICAgICBtZXRob2QgPSAicGVhcnNvbiIsICMgQ29ycmVsYXRpb24gbWV0aG9kIChhbHNvICJzcGVhcm1hbiIgb3IgImtlbmRhbGwiKQ0KICAjICAgICAgICAgICBwY2ggPSAyMSwgICAgICAgICAgICMgcGNoIHN5bWJvbA0KICAjICAgICAgICAgICBsbSA9IFRSVUUsICAgICAgICAgIyBJZiBUUlVFLCBwbG90cyBsaW5lYXIgZml0IHJhdGhlciB0aGFuIHRoZSBMT0VTUyAoc21vb3RoZWQpIGZpdA0KICAjICAgICAgICAgICBjb3IgPSBUUlVFLCAgICAgICAgICMgSWYgVFJVRSwgcmVwb3J0cyBjb3JyZWxhdGlvbnMNCiAgIyAgICAgICAgICAgamlnZ2xlID0gRkFMU0UsICAgICAjIElmIFRSVUUsIGRhdGEgcG9pbnRzIGFyZSBqaXR0ZXJlZA0KICAjICAgICAgICAgICBmYWN0b3IgPSAyLCAgICAgICAgICMgSml0dGVyaW5nIGZhY3Rvcg0KICAjICAgICAgICAgICBoaXN0LmNvbCA9IDQsICAgICAgICMgSGlzdG9ncmFtcyBjb2xvcg0KICAjICAgICAgICAgICBzdGFycyA9IFRSVUUsICAgICAgICMgSWYgVFJVRSwgYWRkcyBzaWduaWZpY2FuY2UgbGV2ZWwgd2l0aCBzdGFycw0KICAjICAgICAgICAgICBjaSA9IFRSVUUpICAgICAgICAgICMgSWYgVFJVRSwgYWRkcyBjb25maWRlbmNlIGludGVydmFscw0KDQoNCmBgYA0KDQpUaGF0J3MgYWxsIQ0KDQohW10oaHR0cHM6Ly9tZWRpYS5naXBoeS5jb20vbWVkaWEveFVPeGY3WGZtcHh1U29kZTFPL2dpcGh5LmdpZikNCg0KDQoNCg0KDQo=